home *** CD-ROM | disk | FTP | other *** search
/ Speccy ClassiX 1998 / Speccy ClassiX 98.iso / amiga_system / the_aminet / dev / gcc / libmat.lha / src / msolve.cc < prev    next >
C/C++ Source or Header  |  1980-01-01  |  1KB  |  75 lines
  1. //                MATRIX LIB
  2. //            TOMMY JOHANSSON 1995
  3.  
  4. #include "matrix.h"
  5. Matrix solve(const Matrix& A,const Matrix& b)
  6. {
  7.     #ifdef DEBUG
  8.         puts("L÷ser ett ekvationssystem.");
  9.     #endif
  10.     
  11.     #ifdef CHECK
  12.         if(A.m!=b.n)
  13.         {
  14.             printf("Felaktiga dimensioner! %d<>%d\n",A.m,b.n);    
  15.             exit(0);
  16.         }
  17.     #endif
  18.  
  19.  
  20.  
  21.     int i,j,o,y;
  22.     float t,p,v;
  23.     Matrix C(A,A.m+1,A.n+1);
  24.     Matrix x(A.m,1);
  25.  
  26.     C.n+=1;
  27.     for(i=1;i<=C.m;i++)
  28.         C.koff[i][C.n]=b.koff[i][1];
  29.  
  30.     
  31.  
  32.     for(o=1;o<=C.m-1;o++)
  33.     {
  34.         for(i=o+1;i<=C.m;i++)
  35.         {
  36.             for(y=o;y<=C.n;y++)
  37.             {
  38.                 if((y==o)&&(C.koff[o][o]!=0)) 
  39.                     v=C.koff[i][o]/C.koff[o][o];
  40.                 else if(C.koff[o][o]==0)
  41.                 {
  42.                     for(i=o+1;i<=A.m;i++)
  43.                     {
  44.                         if(C.koff[i][o]!=0)
  45.                         {
  46.                             for(j=o;j<=A.n;j++)
  47.                                 C.koff[o][j]+=C.koff[i][j];//lΣgg till rad till rad
  48.                             i=A.m;                           // bryt slingan
  49.                         }
  50.                     }
  51.                 }
  52.                 #ifdef DEBUG
  53.                 printf("v=%f o=%d i=%d y=%d\n",v,o,i,y);
  54.                 #endif
  55.                 C.koff[i][y]-=v*C.koff[o][y];
  56.                 
  57.             }
  58.         }
  59.     }
  60.  
  61.     for(i=C.m;i>=1;i--)
  62.     {
  63.         if(C.koff[i][i]==0)
  64.             t=1;
  65.         else
  66.             t=C.koff[i][i];    
  67.     
  68.         p=0;
  69.         for(j=C.m;j>=i+1;j--)
  70.             p+=x.koff[j][1]*C.koff[i][j];
  71.     
  72.         x.koff[i][1]=(C.koff[i][C.n]-p)/t;
  73.     }
  74.     return(x);
  75. }